function [ Out_, Spill_indx, IRF_g, V_decom_g ] = DYFunctReplica(VAR_lag, nMA_lags, Dat_, plot_) 
%%
% Spillover analysis using variance-decomposition, i.e. following 
%    the Diebold-Yilmaz approach. However, in addition to 
%    Cholesky decomposition, I also implement the generalised variance decomposition
%    relying on Peseran-Shin (Economics Letters 1998)
% You can leave here only the Generalized Variance Decomposition
% I m trying to compute here only the Function considering the GVD. 
%  Input:
%  -------------------------------------------------------------------------------
%  VAR_lag  =  number of lags in the VAR
%  nMA_lags =  horizon at which the Spillovers are calculated
%  Dat_     =  the data ( nObs-by-nSeries )
%  plot_    =  -> 1  IRF and VDs are plotted    ->0  no plotting
%
%  Output:
%  -------------------------------------------------------------------------------
%  Out_   contains the spillovers 'from others' and 'to other' for each variable
%         for generalised variance decompositions
%         [ generalised from, generalised to ]
%
%  Spill_indx   is the spillover index  [generalised]
%
%  The IRF are 3D-arrays containing stacked matrices of: 
%       
%       irf( lags, response of variable{j}, shock to variable{i} )
%      

[nObs, nSeries]         = size(Dat_);   %number of series and observations

spec1                   = vgxset('n', nSeries, 'nAR', VAR_lag, 'Constant', true);
[EstSpec, EstStdErrors] = vgxvarx(spec1, Dat_(VAR_lag+1:end,:), [], Dat_(1:VAR_lag,:));  % estimate VAR model
SpecMA                  = vgxma(EstSpec, nMA_lags, 1:nMA_lags );                         % convert VAR to MA representation

SpecMA.MA(2:nMA_lags+1) = SpecMA.MA(1:nMA_lags);
SpecMA.MA{1}            = eye(nSeries);

%
% ...  Impulse Response Functions
%
IRF_g  = zeros(nMA_lags,nSeries,nSeries);
IRF_g1 = zeros(nMA_lags,nSeries,nSeries);

sig_jj = diag(SpecMA.Q);
for ( j=1:nSeries )
    indx_      = zeros(nSeries,1);
    indx_(j,1) = 1;
    for ( k=1:nMA_lags )   % k counts the lag
        IRF_g1(k,:,j) = SpecMA.MA{k}*SpecMA.Q*indx_;  
        IRF_g(k,:,j)  = sig_jj(j,1).^(-0.5).*IRF_g1(k,:,j);                % P-S eqn 10    
    end
end
%
% ...  Variance Decompositions
VD_g1     = zeros(nMA_lags,nSeries,nSeries);
VD_g      = zeros(nMA_lags,nSeries,nSeries);

denom_tmp = zeros(nMA_lags,nSeries);
denom     = zeros(nMA_lags,nSeries);

for ( k=1:nMA_lags )  % diagonal elements of denominator on P-S p.20
    denom_tmp(k,:) = diag(SpecMA.MA{k}*SpecMA.Q*SpecMA.MA{k}');
end
denom = cumsum(denom_tmp,1);

VD_g1   = cumsum(IRF_g.^2);
denom_g = sum(VD_g1,3);
%  dim(IRF) = nHorizon, nSeries, nSeries and is organised as 
%             nHorizon, Response Variables (in each columns), shocked Variables (in each 3d layer) 
%
for (j=1:nSeries)
    VD_g(:,:,j)  = VD_g1(:,:,j)./ denom_g;    % using the definition from Lanne et al
end                                                             

%% Plotting
if (plot_==1)
    figure
    title('IRF - Generalised')
    n=1;
    for ( j=1:nSeries )        % ... x-axis=shocks  and  y-axis  responses
        for ( k=1:nSeries )
            subplot(nSeries,nSeries,n), plot( squeeze(IRF_g(:,j,k)) )
            str_=['Resp ', num2str(j), ' ', ' shk ', num2str(k)];
            title(str_)
            n=n+1;
        end
    end

    % figure
    % title('CUMIRF - Generalised')
    % n=1;
    % for ( j=1:nSeries )        % ... x-axis=shocks  and  y-axis  responses
    %     for ( k=1:nSeries )
    %         subplot(nSeries,nSeries,n), plot( cumsum(squeeze(IRF_g(:,j,k))) )
    %         str_=['CumResp ', num2str(j), ' ', ' shk ', num2str(k)];
    %         title(str_)
    %         n=n+1;
    %     end
    % end
    % 
    % figure
    % title('VD - Cholesky')
    % n=1;
    % for ( j=1:nSeries )
    %     for ( k=1:nSeries )
    %         subplot(nSeries,nSeries,n), plot( 100.*squeeze(VD_oa(:,j,k)) )
    %         str_=['Resp ', num2str(j), ' ', ' shk ', num2str(k)];
    %         title(str_)
    %         n=n+1;
    %     end
    % end

    figure
    title('VD - Generalised')
    n=1;
    for ( j=1:nSeries )
        for ( k=1:nSeries )
            subplot(nSeries,nSeries,n), plot( 100.*squeeze(VD_g(:,j,k)) )
            str_=['Resp ', num2str(j), ' ', ' shk ', num2str(k)];
            title(str_)
            n=n+1;
        end
    end
end
%% Spillover calculation following Diebold-Yilmaz(2009)
%
Spill_g_from = zeros(nSeries,nMA_lags);
Spill_g_to   = zeros(nSeries,nMA_lags);

for (j=1:nMA_lags)
    Spill_g_from(:,j) = sum( squeeze(VD_g(j,:,:)),2)-diag(squeeze(VD_g(j,:,:)));
    Spill_g_to(:,j)   = sum( squeeze(VD_g(j,:,:)),1)'-diag(squeeze(VD_g(j,:,:)));   
end

indx_g = sum(Spill_g_from(:,nMA_lags),1)/(sum(diag(squeeze(VD_g(nMA_lags,:,:))))+sum(Spill_g_from(:,nMA_lags),1) );


Out_       = [ Spill_g_from(:,nMA_lags) Spill_g_to(:,nMA_lags)];
Spill_indx = [indx_g];

V_decom_g  = squeeze(VD_g(nMA_lags,:,:));
